Ground state projection of quantum spin systems in the valence bond basis 
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A Monte Carlo method for quantum spin systems is formulated in the basis of valence bond 
(singlet pair) states. The non-orthogonality of this basis allows for an efficient importance-sampled 
projection of the ground state out of an arbitrary state. The method provides access to resonating 
valence-bond physics, enables a direct improved estimator for the singlet-triplet gap, and extends 
the class of models that can be studied without negative-sign problems. As a demonstration, the 
valence bond distribution in the ground state of the 2D Heisenberg antiferromagnet is calculated. 
Generalizations of the method to fermion systems are also discussed. 
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Quantum spin systems play a prominent role in current 
condensed matter research. An increasing number of ma- 
terials, with a rich variety of lattice geometries and spin 
interactions, realize many different ordered and disor- 
dered quantum states. On the theory side, much progress 
has been made, in particular following the discovery of 
superconductivity in doped antiferromagnetic cuprates, 
but there are still formidable challenges remaining in ex- 
ploring the plethora of possible ground states and excita- 
tions [lj . Gaining deeper understanding of quantum spin 
physics would not only be important in the context of 
particular magnetic systems, but could also give insights 
of broader relevance to correlated quantum matter, e.g., 
concerning quantum phase transitions 0, El • 

Spin models such as the Heisenberg hamiltonian, with 
interactions J^Sj-Sj (with spin Si — 1/2,1, etc.), can be 
studied by a wide range of analytical and numerical meth- 
ods. Since all techniques have limitations, comparisons 
of results obtained in different calculations have proved 
to be crucial. Numerical finite-lattice calculations can 
in principle deliver results free of approximations, but in 
practice available computational methods are restricted 
to certain classes of models. For instance, density matrix 
renormalization £| is essentially limited to one dimen- 
sion and non-approximate quantum Monte Carlo (QMC) 
0, can be used on a large scale only when the in- 
teractions are non-frustrated. Even then, there are still 
often challenges in going to lattice sufficiently large for 
reliable extrapolation to the infinite lattice. Developing 
more general and efficient computational methods, for 
frustrated as well as non- frustrated systems, therefore 
continues to be an important field of research. 

In this Letter, a QMC method for S = 1/2 Heisen- 
berg models is presented which offers several advantages 
relative to state of the art ground state simulations (tak- 
ing T — > in world-line p| or stochastic series expansion 
6] QMC with loop-cluster updates). The method is for- 
mulated in an over-complete and non-orthogonal basis, 
which in the simplest case is the valence-bond (VB) ba- 
sis, in which pairs of spins form singlets 0, 00 E3 Any 
singlet state can be expressed as a sum of VB states, and 



the ground state can be projected out of an arbitrary VB 
state by applying a high power of the hamiltonian H . 
Such a scheme was used by Liang ^lj > who started from 
a good trial wave function | , F) and improved it by sam- 
pling terms of (— H) n \^). However, the projection stage 
apparently did not use importance sampling 12], and 
an extrapolation in n had to be performed. Santoro et 
al. also used the VB basis, employing a Green's function 
method which also does not involve importance sampling 
fliij l . Here it will be shown how the non-orthogonality of 
the VB basis enables a fast importance sampling of the 
terms; no variational state or extrapolations are needed. 
By including triplets, excited states can also be studied, 
and unpaired spins (spinons) can be introduced as well. 
There is thus direct access to degrees of freedom that are 
normally not available with QMC but are of great theo- 
retical interest. The valence bonds and spinons are the 
actors in resonating valence-bond physics Q , which is of- 
ten used as starting point for simplified quantum-dimer 
models and field-theories 0, 01 The method 

should thus facilitate closer contact with modern ana- 
lytical treatments. Moreover, the method extends the 
range of models that can be studied without negative- 
sign problems. 

Projection of a singlet state is here considered first, 
and then the scheme is extended to a triplet. As an illus- 
tration of results that can be obtained, the distribution 
of VB lengths in the ground state of the 2D Heisenberg 
model is presented. Finally, generalizations to a wider 
range of models in other related over-complete bases are 
discussed. 

The expansion in terms of VB states of a singlet ground 
state of N spins (N/2 valence bonds) is written as 

|0> = £ M(aJ,6f) • • • = MSk), (1) 
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where (a^, tff) denotes two spins paired up in a singlet, 



(a,b) = (Uib- Ub)/V2, 



(2) 



i.e., a valence bond, and k labels all bond tilings of the 
lattice (allowing arbitrary bond lengths). The notation 
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FIG. 1: Action of a bond operator on two VB states. 
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FIG. 2: Singlet convention for a system with sublattices A and 
B. The arrows indicate the order of the spins in the singlets, 
e.g., c^b means (c, b) = (felt - |cTb)/v / 2- 



\Sk) has been introduced for convenience. The expansion 
can always be made positive definite; any negative fk can 
be made positive by switching the indices of one singlet. 

Since the VB basis is over-complete, the expansion co- 
efficients are not unique in general. However, the expan- 
sion of any state I'F) written in the VB basis of course 
has a unique expansion in energy eigenstates |n); 

k n 

Therefore, acting on this state with a high power (n — > 
oo) of the hamiltonian projects out the ground state: 

[-(H-C)] n \V) =Y,9n,k\Sk) ^c \E -C\ n \0). (4) 

k 

A constant C has been subtracted in order to render the 
magnitude of the lowest eigenvalue larger than the high- 
est one. QMC methods based on in the standard 
basis of eigenstates of all Sf are commonly used 0, 0] , 
although for bipartite systems they tend to be less ef- 
ficient than low-temperature simulations with advanced 
finite- T methods la la- However, for frustrated systems 
[l7j and t-J models |l8j . where there are sign problems 
(non-positive-definite g Ul k and analogous mixed signs in 
other methods) projector methods are superior if a good 
trial state can be used. 

The first observation underlying the projector method 
in the VB basis is that the application of a Heisenberg 
interaction operator on a VB state leads to a very simple 
rearrangement of valence bonds 0, |n|. Consider the 
Heisenberg hamiltonian, on any lattice, written as 

H = - 2J JijHij, H i:j = -(Si ■ Sj ~ ~). (5) 

Acting with H a b on a VB state in which sites a and b 
belong to the same valence bond gives an eigenvalue of 
unity; H a b\..(a, &)..) = |.. (a, &)..). Acting on sites belong- 
ing to different valence bonds gives a new basis state; 

H bc \..(a,b)..(c,d)..) =i|..(M)..(e, &)..). (6) 

This bond flip is illustrated in Fig. LI Sites a and d 
are completely arbitrary and the sign is always positive 
when the indices are in the order indicated. This im- 
plies that a positive-definite representation of the projec- 
tion {—H) n \9) can be achieved for a bipartite lattice, by 



defining (a, b) so that a is always on sublattice A and b is 
on sublattice B. This is illustrated using arrows on the 
bonds in Fig. |3 The convention also implies a positive- 
definite expansion [Io| . 

For a non-frustrated interaction, one can show that the 
projected ground state (@J contains only bonds connect- 
ing spins on different sublattices (bipartite bonds) . Con- 
sider the VB configuration shown to the left in Fig.QJb). 
With sites a,c in sublattice A and b,d in B, both the 
bonds are non-bipartite. When the operator has acted, 
the new bonds are bipartite. With a bipartite interaction 
one cannot accomplish the reverse process (note that the 
hamiltonian is not manifestly hermitean in the VB ba- 
sis) and thus, if the trial state |W) contains non-bipartite 
bonds, they will vanish after H has acted a number of 
times and cannot reappear. The ground state hence must 
contain only bipartite bonds. 

For a frustrated interaction, non-bipartite bonds are 
generated and the bond flip © can lead to a minus sign 
La. it may also in practice not be possible to find a 
singlet convention which renders all the expansion co- 
efficients in {Q positive. Nevertheless, the projection 
scheme in principle works also for frustrated systems, as 
long as any negative signs are taken into account in the 
standard way |2Qj. 

With the hamiltonian C = can be used in (@J. To 
expand H n , an index sequence P n — [<xi, £>i], . . . , [a n , b n ] 
is used to refer to an operator product H ap b p . One of 
the VB basis states can be chosen as the trial state; |W) = 
\So) HD< The projected state is then 

n 

(-H) n \S ) = ^2l[Ja pbp H apbp \So) 

n 

= Ell^l^n)}, (7) 

Pn P=l 

where \S(P n )) denotes the (normalized) state obtained 
when the operators have acted on |So). The factors 
w a i, are \ J a b or J a j, for operations with H a b that cause 
bond flips and are diagonal, respectively, in the course of 
propagating the state from |S ) to \S(P n )). The weights 
W(P n ) = Y[ p w a p b p can be identified with the expansion 
coefficients g n ^ in J3J. Note that no operator H a b can de- 
stroy the state in 0, i.e., all W(P n ) ^ 0. This "softness" 
of the basis will be taken advantage of in constructing an 
efficient importance-sampling scheme that would not be 
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FIG. 3: Two VB states in two dimensions and their overlap 
in terms of loops formed by superimposing the two bond con- 
figurations. In this case there are N v = 8 valence bonds and 
Ni = 3 loops, and (S a \S ) = 2 N '- N " = 1/32 0. 

possible in an orthogonal ("hard") basis, where there are 
enormous constraints on the operator products. 

To calculate the ground state energy, the overlap with 
an arbitrary reference state \R) can be taken; 

(R\H\0) Ep„ W(P n )(R\H\S(P n )) 
° (R\0) Ep„ W(P n )(R\S(P n )) ■ [ > 

One can always choose a state with equal overlap with 
all basis states, and hence all (R\S) overlaps cancel. If 
P n is sampled with probability cx |W(P„)|, the energy is 

Eo = -j-- Jij (s(rHj + |(1 - rii^qij)) , (9) 

where riij = 1 (0) if there is (is not) a bond connecting 
sites i and j in \S(P n )). For a frustrated system, s = ±1 
is the product of phase factors ±1 arising when propa- 
gating | So) t° \S(P n )) with the string P n , and q^ = ±1 
arises when Hij is applied once more. For a bipartite 
system, s = q io = 1 and thus Eq = — § Ji.j( n ij + !)■ 
An expectation value (A) = (0|A|0)/(0|0) of an arbi- 
trary operator can be written in terms of two projected 
states, obtained from the same trial state \ Sq) propagated 
with two different operator strings P n and Q„; 

, , _ Z Pn Z Qn W(P n )W(Q n )(S(Q n )\A\S(P n )) 

{ } E P ,XQ n W ( P ^ W (Qn)(S(Q n )\S(P n )) ■ [ Uj 

The weight function to be used in importance sampling 
is thus W(P n )W{Q n )(S(Q n )\S(P n )), and the operator 
estimator is (S(Q n )\A\S(P n )) / (S(Q n )\S(P n )) . 

For a bipartite system, the overlap of two VB states is 
determined by the loops formed when the bonds are su- 
perimposed [22j, as illustrated in Fig. [3| (with frustration, 
a sign has to be determined as well). Matrix elements 
(S Q |S, • Sj\Sp) are also easily obtained from these loops; 
(S a \S t ■ Sj|S^)/(5 a |S» = ±3/4 if sites i and j belong 
to the same loop (+ and — for i, j on the same and dif- 
ferent sublattices, respectively, in the case of a bipartite 
lattice), and otherwise. 

A remarkable aspect of the VB basis is that Eqs. © 
and (|1C)|> can be efficiently sampled in an almost triv- 
ial way, in steps where a few (r) of the operators in the 
product P n are changed at random. Naively one might 



expect that the acceptance rate should become very low 
for large expansion order n, but this turns out not to 
be the case. With r = 4, the acceptance rate in the 
case of the 2D Heisenberg model is m 40%, almost in- 
dependently of n and the lattice size. The new weight 
can be computed by performing the full propagation of 
the state \Sq) with the updated product(s) in (J7J) [and 
calculating the new overlap in the case of <|10jl]. Recalcu- 
lating the full weight, instead of just a ratio, may seem 
like an inefficient proposition. However, if n has to be in- 
creased with the systems size as N a in order to converge 
to the ground state, N 2a operations [N 1+a if a < 1 in 
the case of <|10[l ] are needed to update the full operator 
sequence (attempting n/r updates of r operators is de- 
fined as one sweep; several measurements are carried out 
during each sweep). In T — * calculations with finite-T 
methods dEI the scaling is N 1+a ' if T cx N~ a ' . Hence, 
if a, a! rs 1 the scaling is very similar. The gap to the 
lowest singlet excitation dictates a and hence in many 
cases a < 1 suffices. An even faster sampling could likely 
be achieved by using a linked operator list 0; such an 
improvement will be left for future work. 

In order to study a triplet state, consider a triplet 
bond; [a, b} = (taifc + latb)/V^- The eigenvalue of H ab 
operating on [a, b]o is 0. If Hb c is applied to [a, 6]q(c, d), 
the reconfiguration of the bonds is exactly as in JBJ; the 
new state is [a, d]o(c, b)/2. Hence, if there is no diagonal 
operation on the triplet, the triplet bond (if there is only 
one) behaves exactly as a singlet and the only change in 
the scheme is in the operator estimators. This enables 
an improved estimator for, e.g., the singlet-triplet gap: 
Carrying out the simulation with only singlets, one of 
the bonds can be flagged as a triplet at the measurement 
stage. The E\ estimator can be averaged over all N/2 ini- 
tial triplet choices, with contributions coming only from 
surviving configurations, i.e., those for which there are 
no diagonal operations on the triplet (the survival ratio 
depends on n). This does not change the scaling N 2a 
of the simulation and can vastly improve the estimate of 
the gap compared to E\ — Eq obtained from two inde- 
pendent simulations (the improvement is mainly due to 
partial cancelation of correlated statistical errors in Eq 
and Ei). For example, for the 2D Heisenberg model with 
N = 64x 64, a projection with n = 15N and 10 6 updat- 
ing sweeps gave Eq/N = 0.669449(2) and the finite-size 
gap Ei — Eq = 0.0041(2), corresponding to an accuracy 
gain of 60 times for the gap, or a CPU-time reduction 
of 7000. The energy agrees with that obtained using the 
SSE method @; E„/N = 0.669450(1), confirming the 
unbiased nature of both calculations. 

It is important to verify that the method works for 
frustrated interactions as well, although the basic formu- 
lation discussed here 0] is not likely to be practically 
useful for large frustrated lattices. The method should, 
however, be applicable to models with local sign prob- 
lems, e.g., frustrated impurities (in large host lattices). 
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FIG. 4: The bond-length probability for bonds along the line 
(x, 0) in the 2D Heisenberg ground state wave function. 



Checks against exact diagonalization results confirm that 
the scheme indeed works. For a 4 x 4 system with nearest- 
and next-nearest-neighbor interactions J\ and J2, at 
J 2 /J 1 = 0.1 a ground state energy E a /N = 0.65986(4) 
was obtained using n = 3iV and \Sq) a columnar dimcr 
state (5 x 10 9 updating sweeps), which matches well the 
exact Eo/N = 0.659817. The average sign in this case is 
(s) « 0.074. 

As an application of the method, the VB length distri- 
bution in the ground state of the 2D Heisenberg model is 
presented next. Liang et al. 10] studied variational wave 
functions with VB state amplitudes fk — Jli M a ?> ) m 
Eq. JQ). For h cx l/r p , where r is the bond length, they 
concluded that there is long-range Neel order for p < 5. 
The best variational energy was obtained with p = 4, 
but the dependence on p for 2 < p < 5 was quite weak. 
The bond amplitude h(x, y) does not correspond exactly 
to the probability P(x,y) cx Y^ fe fkn^y [where n^. y is the 
number of length-(a;, y) bonds in VB state fc], but Monte 
Carlo simulations of the type used in Ref. [li| confirm 
that if h(r) cx l/r p then also P(r) cx l/r p . Note also 
that P(x, y) is not a ground state expectation value but 
a property of the wave function coefficients [but the ex- 
pectation value (Oln^ylO) turns out to be almost identical 
to P(x,y)]. A potential worry is that since the VB ba- 
sis is over-complete, the bond distribution is not unique. 
However, the way the projection is done corresponds to 
a uniform averaging over all possible VB ground state 
wave functions; P(x,y) defined this way clearly has a 
well-defined meaning. 

Calculations were carried out on periodic Lx L lattices 
with L = 64 and 128, with n up to 15iV and 20iV, respec- 
tively (convergence was checked). The results shown in 
Fig. 0] suggest that P(r) ~ 1/r 3 (there is no notable an- 
gular dependence) . It would be interesting to understand 
this result from an analytical starting point, and also to 
study the probability distribution for a quantum-critical 
system, e.g., the Heisenberg bilayer [23j. 

The projector scheme discussed here opens up a range 



of interesting and promising avenues to be explored. The 
VB (+ triplets) basis is formed out of the 2-site eigen- 
states of Si ■ Sj and hence is particularly suitable for 
Heisenberg models. The method can also be extended, 
without sign problems, to higher-order (non- frustrating) 
interactions of the form — (Sj • Sj — j)(Sk • Sj — j). This 
interaction has a sign problem in the z-basis, and hence 
the present method solves a class of sign problems. Al- 
though there are sign problems for frustrated systems in 
general, the VB basis opens opportunities to explore can- 
celation schemes based on over-completeness [l9|. Good 
sign-problem-free approximations could also perhaps be 
developed. 

It is possible to generalize the VB basis to other Hilbert 
spaces, with different types of bonds corresponding to 
eigenstates of the Hamiltonian on two sites (or even > 2 
sites, although the complexity of the approach then in- 
creases considerably). Such schemes for t-J and Hubbard 
models will be investigated. Although there will clearly 
be sign problems for fermions, a generalized bond-state 
basis also offers opportunities for new variational wave 
functions, which could be further refined with the pro- 
jector method. Access to the VB (and generalized) de- 
grees of freedom also enables construction of interesting 
hamiltonians acting on bonds. Such studies could clarify 
the relationships between quantum dimer |l4j and spin 
models. Studying the properties (such as the length- 
distribution) of a triplet bond in the "singlet soup" of 
a gapped or critical system gives information pertain- 
ing to the nature of the spinon bound state (magnon) or 
spinon doconfinement. This should be very useful, e.g., 
in studies of deconfined quantum-criticality [si Il5j|. 
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